/*
 * unconstrained.c
 *
 *  Created on: Jun 19, 2012
 *      Author: wlin15
 */

#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include "qrmcp.h"
#include <gsl/gsl_permutation.h>
#include <math.h>
     int
     main (void)
     {

    	 size_t i,j;
/*    	 const size_t M = 6,N=6;
    	 gsl_matrix *A = gsl_matrix_alloc (M, N);

     	 gsl_vector * y = gsl_vector_alloc (M);

    	 //Initialize A
    	 for (i = 0; i < M; i++)
    		 for (j = 0; j < N; j++)
    			 gsl_matrix_set (A, i, j, i*i*i+j*j*j);

    	 //Initialize y
    	 for (i = 0; i < M; i++)
    	 {
    		 gsl_vector_set (y, i,  i*i);
    	 }*/

    	 size_t m;
    	 m=6;
    	 gsl_permutation * piv = gsl_permutation_alloc (m);
    	 gsl_matrix *A=gsl_matrix_alloc(m,m);
    	 gsl_vector *y=gsl_vector_alloc(m);
         FILE *fA=fopen("A.dat","rb");

         FILE *fy=fopen("y.dat","rb");

       //gsl_matrix_fread(fA,A);
         gsl_matrix_fscanf(fA,A);
         gsl_vector_fscanf(fy,y);
         fclose(fA);
         fclose(fy);
    	 //Initialize permutation vector
		 gsl_permutation_init (piv);


        //print out the matrix A
		 printf("The original matrix A is:\n");
		 for (i = 0; i < m; i++)
		 {
			 for (j = 0; j < m; j++)
			 {
				 printf("%e",gsl_matrix_get (A, i, j));

			 }
			 printf("\n");
		 }


/*		 printf("\n");

		 //print out the vector y
		 printf("The original y is\n");
		 gsl_vector_fprintf(stdout,y,"%e");

         printf("\n");
		 //print out the original permutation vector
		 printf("The original permutation vector piv is\n");
		 gsl_permutation_fprintf (stdout, piv, " %12zu");*/


       qrmcp (A, y, piv);
       gsl_matrix_free (A);
       gsl_vector_free(y);
       gsl_permutation_free(piv);

       return 0;
     }
